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Abstract —Quantitative analysis of channel networks plays 
an important role in river studies. To provide a quantitative 
representation of channel networks, we propose a new method 
that extracts channels from remotely sensed Images and estimates 
their widths. Onr fnlly automated method is based on a recently 
proposed Multiscale Singularity Index that responds strongly 
to cnrvillnear structures hut weakly to edges. The algorithm 
produces a channel map, nsing a single image where water 
and non-water pixels have contrast, snch as a Landsat near- 
infrared hand image or a water index defined on mnltiple 
hands. The proposed method provides a robust alternative 
to the procednres that are used in remote sensing of fluvial 
geomorphology and makes classification and analysis of channel 
networks easier. The sonrce code of the algorithm is available at: 
http://live.ece.ntexas.edu/researcti/cne/, 

Index Terms —Channel network extraction, river width, deltas, 
remote sensing, image processing. 

I. Introduction 

A method for completely automatic extraction of channel 
networks from satellite imagery could greatly facilitate 
the monitoring of water resources by eliminating the labo¬ 
rious process of manual inspection. Such a method could 
be used for creating quantitative representations of channel 
networks, which would be useful in a wide variety of studies. 
The automatic extraction of channel networks is particularly 
challenging in coastal areas due to low topographic gradients, 
the presence of features such as sediment plumes, and the 
wide range of scales over which channel features are present. 
A robust channel extraction method would ease monitoring 
coastal areas and analyzing deltaic response to anthropogenic 
and natural forcings over large spatial areas and long temporal 
intervals. 

Several approaches have been suggested to extract curvi¬ 
linear structures from remotely sensed images, many of them 
focusing on road network extraction Q-©- Based on the road 
network extraction work 0. 0. a method has been proposed 
to detect rivers as linear structures, imposing constraints on 
river length, curvature, and confluences for connectivity Q . A 
software tool, RivWidth Q, has been proposed for calculating 
river centerlines and widths. RivWidth (vO.4) requires the 
availability of a previously defined binary mask that indicates 
water and non-water pixels. Although such a mask could 
be extracted from remotely sensed images by thresholding 
and using shape-correcting operations, manual cleaning of the 
mask would often be necessary to separate the true water mask 
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from spurious responses ©. Our new method can estimate 
channel centerline, width, and orientation, and create a map 
of a channel network in a purely automatic manner, using only 
remotely sensed data. To the best of our knowledge, this is the 
first fully automatic approach that provides these outputs. 

II. Modified Multiscale Singularity Index 
The Multiscale Singularity Index lTg,|TT) is a useful method 
for detecting singular curvilinear structures over multiple 
scales. The algorithm is useful for locating channels in satellite 
images. However, presence of channels over a wide range of 
scales creates some artifacts in the Singularity Index response. 
Our work modifies and extends the Multiscale Singularity In¬ 
dex to address the multi-scale nature of channel networks. The 
Multiscale Singularity Index algorithm and our modifications 
are explained briefly in the following sections. 


A. The Multiscale Singularity Index 

At each pixel, the Multiscale Singularity Index algorithm first 
estimates the direction 9 orthogonal to the curvilinear mass, 
using second order derivatives of an input image along evenly 
spaced directions. Then, it computes the Singularity Index at 
each scale as; 


{'ipf){x,y,a) = 


\fo,e,crix,y)f2,9,crix,y)\ 


( 1 ) 


1 + \fl,9,aaix,y)\ 

where fo,e,a{x,y), fi^e,a{x,y), and f 2 , 9 ,a{x,y) are the 
zero, first, and second order Gaussian derivatives at scale a 
and along the direction 9{x,y). In the denominator, a is a 
constant with a recommended value 1.7754. This value results 
in maximum attenuation of the side lobe response of the index 
|10|. The index ('(/i/)(x,y,cr) is computed over N scales 
^ for n = 1,2,...,V. The window sizes for 
the image Alters are determined to be [6cr]. Since a channel 
of width larger than the image dimensions cannot be detected 
by the algorithm, an upper bound for the number of scales A 
can be determined by having the Alter dimensions smaller than 
the image dimensions M x M so that 6ai\/2^^ < M: 


N = 


21og|^ 

log 2 


+ 1 


( 2 ) 


After computing the Singularity Index at each scale, the 
algorithm finds the maximum response across all scales at 
each pixel location. 

The Singularity Index retains polarity, which is useful for 
discriminating between channel and island response, since 
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channels and islands have opposite polarity. By discarding the 
negative polarity, we remove the island response. 


B. Debiasing Input Images 

An input image should be debiased before computing the 
Singularity Index in order to achieve invariance to local DC 
offset. The Multiscale Singularity Index algorithm debiases 
an input image by subtracting a large Gaussian filter from 
the original image, which essentially performs local mean 
subtraction. This approach works well for a small range of 
scales. For a large range of scales, however, a large Gaussian 
filter fails to debias finer scales and results in a loss of detail 
at fine scales (Fig. lai. To address this problem, we debias 
the input image at every scale. Instead of using one large 
Gaussian filter over all scales, our modified version of the 
Multiscale Singularity Index uses a Gaussian filter with a 
standard deviation of cr„ at each scale as: 


Ia= I-Qa*I (3) 

where I^r is the debiased image at scale cr, I is the input image, 
and is the Gaussian filter. 


C. Width Estimation 

The channel width is estimated by interpolating between the 
scale that has the highest Singularity Index response '0 and its 
neighbor scales as follows: 

+ 1 

w{x,y) = -— (4) 

E {i^f){x,y,(Tm+z) 

i=-l 

where m = argmax„('0/)(a:, y, an) at spatial coordinate x, y 
and k is a scalar variable that scales the output. 


D. Adaptive Smoothing 

The Multiscale Singularity Index creates ripples near the banks 
of wide rivers when a large range of scales is processed (Fig 


Ici. The ripples occur after finding the maximum response 


over all scales at each spatial coordinate. To attenuate the 
ripples, we employ an adaptive smoothing algorithm that 
adjusts the strength of smoothing based on the estimated scale 
for each pixel, so that coarse scales can be smoothed more than 
fine scales. To implement an adaptive smoothing algorithm in 
a computationally efficient way, we first compute an integral 
image over the Singularity Index response, which enables fast 
computation of summations over regions of arbitrary size. 
Then, we smooth the response using a box filter with a variable 
window size that is determined by the estimated scale. Since 
the integral image is computed only once, the algorithm only 
needs to perform two additions and one subtraction per pixel. 
We apply the adaptive box filter iteratively to approximate a 
Gaussian filter. 

Typical results delivered by the Multiscale Singularity Index 
and our modified version are compared in Fig. [T] 




(C) (d) 


Fig. 1. Comparison of original and modified Multiscale Singularity Index 
response. Full sized images: original response (a) and modified response (b). 
Zoomed images: original response (c) and modified response (d). 


E. Centerline Extraction 

To determine the channel centerlines, the maximum response 
across all scales is computed at each coordinate, and the 
orientation value at the maximum-response scale is taken 
to be the dominant centerline direction. A process of non¬ 
maxima suppression is applied along the dominant direction 
as explained in ng- Then, a threshold level T is determined 
on the non-maxima suppressed (NMS) image using Otsu’s 
method HI- To preserve channel connectivity, a hysteresis 
threshold is applied to binarize the NMS image, as follows: 

1) Set pixels above an error threshold e to one and the rest 
to zero 
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2) Find the connected components in the image 

3) Find and remove the connected components that do not 
have at least one pixel above the threshold T. 

The error threshold e is chosen empirically as 0.1 x T. Ex¬ 
tracted centerlines for an example input image are illustrated 
in Fig. 1^ The hgure is inverted for better visualization. 

F. Creating a Map of Channels 

To show the computed channel width at each spatial coordinate 
along a centerline, a map of channels is created by re-growing 
the channels. The channels are re-grown by drawing a line of 
length w{x,y) and orientation 6{x,y) at each spatial location 
{x,y). The algorithm estimates the width of any water body 
of width smaller than the largest scale. Therefore, the resultant 
map includes ponds and other small water bodies as well as 
channels. Results of channel map creation are presented in 
Section [111] 

III. Experiments and Results 

We tested our model on three different regions having 
different characteristics: Mississippi river near Memphis, 
TN (II), Mississippi, Missouri, and Illinois rivers near St. 
Louis, MO (12), and a portion of the Ganges-Brahmaputra- 
Jamuna Delta (13). We used Landsat-8 images, downloaded 
at http://earthexplorer.usgs.gov/, to create the input images for 
our algorithm. The algorithm requires input images to have 
a contrast between water and non-water pixels. An example 
input could be a near-infrared image or a water index that 
uses multiple bands. In our experiments, the input images were 
created using the Modified Normalized Difference Water Index 
(MNDWI) p3) , which is an effective way to extract water 
information from remote sensing imagery. 

The ground truth for II and 12 are obtained by aligning, 
cropping, and rasterizing the river data from the National Fly- 
drography Dataset (http://nhd.usgs.gov/). Eor 13, the Ganges- 
Brahmaputra-Jamuna (GBJ) Delta network extracted by |9| 
is used as the ground truth. The extraction performed by 191 
included manual cleaning and a comparison to Google Earth 
imagery. 


We fixed the minimum scale cti to its default value 
1.5 pixels, which is the smallest width for a channel to be 
captured by the algorithm. The number of scales is determined 
automatically using the upper bound that is described earlier. 
To reduce the computation time, a smaller number of scales 
could also be chosen if all channels of interest are known to 
be smaller than a certain width. In the experiments reported 
here, we set the number of scales N to 16. 

The regrown channel maps, showing the estimated location, 
width, and orientation of the channels, are compared with the 
ground truth and input images in Eig. The ground truth im¬ 
ages did not include non-channel water bodies. Therefore, we 
removed the non-channel water bodies also from the regrown 
channel maps, by discarding the connected components that 
constitute less than 0.1% of the maps. Given the ground truth 
images, the accuracies ({TP+TN)/{TP + TN + FP + FN)) 
of the re-grown channel network images were found to be 
96.77%, 97.86%, and 91.13% for II, 12, and 13, respectively. 

IV. Conclusion AND Euture Work 

We have described an automatic channel network extraction 
algorithm that inputs remotely sensed images and produces 
maps of estimated channel centerline, width, and orientation. 
We modihed a Multiscale Singularity Index to extract a 
network of channels over a wide range of scales. The algorithm 
works automatically without any user intervention. 

Our method can be used to analyze channel networks in 
different environments and over time to capture the effect of 
environmental forcing and natural and anthropogenic change 
on the network structure. One of our future research directions 
is to analyze deltaic response to anthropogenic and natural 
forcings in coastal areas. We also plan to extend our work 
towards automatically creating topological maps, which will 
provide graph representations of channel networks. 
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